clear all
set more off
set varabbrev off
clear matrix
clear mata
set maxvar 20000
cd "V:\Projects\Enum_Exp_DD\Replication_files"									//Update the path to match your local directory
global root "V:\Projects\Enum_Exp_DD\Replication_files"							//Update the path to match your local directory
capture mkdir "$root\Tables"
global tables "$root\Tables"
capture mkdir "$root\Figures"
global figures "$root\Figures"
version 17


/*==============================================================================
					Appendix 4:Respondent Balance
==============================================================================*/

use "$root\Data\hh_dsy_cr_rec.dta", clear


egen asset_count = rowtotal(sec2_q12_1 sec2_q12_2 sec2_q12_3 sec2_q12_4 sec2_q12_5 sec2_q12_6 sec2_q12_7 sec2_q12_8 sec2_q12_9 sec2_q12_10 sec2_q12_11 sec2_q12_12 sec2_q12_13 sec2_q12_14 sec2_q12_96)
lab var asset_count "count of hh assets"



la var sec3_q1 "Household Size (members)"
la var sec2_q16 "Household has access to electricity"
la var m_uas_inc "Mother's Income"
la var sec1_q7 "Father living in same HH"
la var scst_caste "Whether Backward Caste (SC/ST)"
la var hfc_health "Annual Health Expenditure"
lab var asset_count "Count of Household Assets"
la var prop_educated " Proportion of educated members"



	
// Export to Excel
ds sec3_q1 sec2_q16 m_uas_inc sec1_q7 scst_caste hfc_health asset_count prop_educated

iebaltab `r(varlist)', ///
    grpvar(group_code) ///
    savexlsx("$tables/Appendix_4.xlsx") ///
    rowvarlabels ///
    feqtest ///
    replace ///
    total ///
    order(1 2 3) ///
    addnote("Note: Standard errors are in parentheses. The F-test for joint orthogonality accounts for experimental block fixed effects.")



ds sec3_q1 sec2_q16 m_uas_inc sec1_q7 scst_caste hfc_health asset_count prop_educated

// Export to LaTeX
ds sec3_q1 sec2_q16 m_uas_inc sec1_q7 scst_caste hfc_health asset_count prop_educated

iebaltab `r(varlist)', ///
    grpvar(group_code) ///
    savetex("$tables/Appendix_4.tex") ///
    rowvarlabels ///
    feqtest ///
    replace ///
    total ///
    order(1 2 3) ///
    addnote("Note: Standard errors are in parentheses. The F-test for joint orthogonality accounts for experimental block fixed effects.")




/*==============================================================================
					Appendix 9 (Figure A1-Figure A5)
==============================================================================*/


/*==============================================================================
			Figure A1: (Coefficient Plot of Impact on Data Quality
==============================================================================*/



use "$root\Data\hh_dsy_cr_rec.dta", clear

//Variable Preparation

// Create dummy variables for categorical variables 
// Religion
tab sec2_q1, gen(religion_)

// Caste - particularly important for social stratification
tab sec2_q2, gen(caste_)

// Household ownership type
tab sec2_q6, gen(house_own_type_)

// Direct Benefit Transfer (DBT) receipt
gen dbt_receipt = (sec2_q4==1)
label var dbt_receipt "Household receives DBT"

// Water source - improved vs. unimproved
gen improved_water = inlist(sec2_q18, 9, 10, 13, 14)
label var improved_water "Access to improved water source"

// Sanitation
gen toilet_ownership = (sec2_q20==1)
label var toilet_ownership "Household owns toilet"

gen improved_sanitation = inlist(sec2_q21, 2, 3, 4, 5) if sec2_q20 == 1
label var improved_sanitation "Access to improved sanitation"

// Electricity access
gen electricity = (sec2_q16==1)
label var electricity "Has electricity access"

// Internet access
gen internet = (sec2_q17==1)
label var internet "Has internet access"

// Cooking fuel
gen clean_cooking = 0
replace clean_cooking = 1 if sec2_q24_1 == 1 | sec2_q24_2 == 1 | sec2_q24_3 == 1
label var clean_cooking "Uses clean cooking fuel"

// Child mortality indicator
gen child_mortality = (sec2_q23==1)
label var child_mortality "Child mortality in past 5 years"

// Working mother
gen mother_working = (sec6_q3==1)
label var mother_working "Mother worked in last 365 days"

// Create asset index/score
egen asset_count = rowtotal(sec2_q12_1-sec2_q12_14 sec2_q12_96)
gen asset_proportion = asset_count / 15
label var asset_proportion "Proportion of possible assets"
egen asset_zscore = std(asset_proportion)
label var asset_zscore "Standardized asset score"

// Maximum education in household
egen max_educ_hh = rowmax(sec3_q11_*)
label var max_educ_hh "Maximum education level in household"

// Household size
gen hhsize = sec3_q1
label var hhsize "Household size"

// Number of rooms (proxy for house size)
gen num_rooms = sec2_q7
label var num_rooms "Number of dwelling rooms"

// Property value
gen property_value = sec2_q8
label var property_value "Total property value"

// DBT amount
gen dbt_amount = sec2_q5
label var dbt_amount "DBT amount received"

// Health expenditure
egen total_health_exp = rowtotal(sec3_q18_*), missing
label var total_health_exp "Total health expenditure (12 months)"

// Education expenditure
egen total_edu_exp = rowtotal(sec3_q17_*), missing
label var total_edu_exp "Total education expenditure (12 months)"

// Number of siblings under 12
gen siblings_under12 = sec7_q2
label var siblings_under12 "Number of siblings under 12"

label var asset_zscore "Asset index"
label var improved_sanitation "Improved sanitation facility"
label var improved_water "Improved water source"
label var max_educ_hh "Maximum education in household"


// HOUSEHOLD CONTROLS - DUMMY
global hh_d_contr "religion_* caste_* house_own_type_* dbt_receipt electricity internet improved_water toilet_ownership clean_cooking child_mortality mother_working"

// HOUSEHOLD CONTROLS - CONTINUOUS
global hh_c_contr "hhsize max_educ_hh num_rooms property_value total_health_exp total_edu_exp asset_zscore siblings_under12"



tab c_strata, gen(c_strata_)
// STRATA CONTROLS 
global str_contr "c_strata_*"
global str_contr "i.c_strata"




**FINAL CONTROL VARIABLE SET FOR PDS LASSO, we  need to replace the missing value as recommended by the David Mckinzie 

global miss_hh_d_contr ""
global dmiss_hh_d_contr ""

foreach var of varlist $hh_d_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_d_contr "$miss_hh_d_contr miss_`var'"
    global dmiss_hh_d_contr "$dmiss_hh_d_contr dmiss_`var'"
}

global miss_hh_c_contr ""
global dmiss_hh_c_contr ""

foreach var of varlist $hh_c_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_c_contr "$miss_hh_c_contr miss_`var'"
    global dmiss_hh_c_contr "$dmiss_hh_c_contr dmiss_`var'"
}


global controls "$str_contr $miss_hh_d_contr $dmiss_hh_d_contr $miss_hh_c_contr $dmiss_hh_c_contr"




zscore outliers vskips invalid_skips_c inconsistency

eststo clear

// Run the same regressions as before to get the estimates - USING Z-VARIABLES
foreach var of varlist whether_outlier z_outliers z_vskips z_invalid_skips_c z_inconsistency {
    areg `var' t1 t2 i.enum_id, a(c_strata) vce(cluster vill_id)   
    eststo ols_`var'
}

foreach var of varlist whether_outlier z_outliers z_vskips z_invalid_skips_c z_inconsistency {
    pdslasso `var' t1 t2 ($controls), cluster(vill_id) partial($str_contr)
    eststo pds_`var'
}



// 1. Whether Data Quality Flag 
#delimit ;
coefplot 
    (ols_whether_outlier, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
        msymbol(circle) msize(medium) label("OLS"))
    (pds_whether_outlier, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
        msymbol(diamond) msize(medium) label("PDS-Lasso")),
    
    vertical
    xtitle("")
    xscale(range(0.5 2.5))
    xlabel(
        1 "T₁" 
        2 "T₂", 
        labsize(medium) angle(0)
    )
    
    yline(0, lpattern(dash) lcolor(gs10))
    yscale(range(-0.6 0))   
    ylabel(-0.6(0.3)0, format(%9.2f) labsize(small))
    ytitle("Treatment Effect", size(medium))
    title("a) Whether Outliers (z-score)", size(medsmall) color(black))
    
    legend(off)
    
    mlabel format(%9.2f)
    mlabsize(medium)          
    mlabposition(3 9 3 9)     
    mlabgap(*1.5)             
    mlabcolor(black)          
    
    graphregion(color(white) margin(small))
    bgcolor(white)
    plotregion(color(white) margin(medium))
    name(plot_direct_1, replace);
#delimit cr

// 2. Data Quality Flags (Outliers) - USING Z-VARIABLES
#delimit ;
coefplot 
    (ols_z_outliers, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
        msymbol(circle) msize(medium) label("OLS"))
    (pds_z_outliers, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
        msymbol(diamond) msize(medium) label("PDS-Lasso")),
    
    vertical
    xtitle("")
    xscale(range(0.5 2.5))
    xlabel(
        1 "T₁" 
        2 "T₂", 
        labsize(medium) angle(0)
    )
    
    yline(0, lpattern(dash) lcolor(gs10))
    yscale(range(-1.2 0.4))
    ylabel(-1.2(0.4)0.4, format(%9.2f) labsize(small))
    ytitle("Treatment Effect", size(medium))
    title("b) Outliers (z-score)", size(medsmall) color(black))
    
    legend(off)
    
    // Place coefficient values to the sides of the markers
    mlabel format(%9.2f)
    mlabsize(medium)         
    mlabposition(3 9 3 9)    
    mlabgap(*1.5)            
    mlabcolor(black)         
    
    graphregion(color(white) margin(small))
    bgcolor(white)
    plotregion(color(white) margin(medium))
    name(plot_direct_2, replace);
#delimit cr

// 3. Valid Skips
#delimit ;
coefplot 
    (ols_z_vskips, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
        msymbol(circle) msize(medium) label("OLS"))
    (pds_z_vskips, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
        msymbol(diamond) msize(medium) label("PDS-Lasso")),
    
    vertical
    xtitle("")
    xscale(range(0.5 2.5))
    xlabel(
        1 "T₁" 
        2 "T₂", 
        labsize(medium) angle(0)
    )
    
    yline(0, lpattern(dash) lcolor(gs10))
    yscale(range(-0.9 0.3))
    ylabel(-0.9(0.3)0.3, format(%9.2f) labsize(small))
    ytitle("Treatment Effect", size(medium))
    title("c) Valid Skips (z-score)", size(medsmall) color(black))
    
    legend(off)
    
    // Place coefficient values to the sides of the markers
    mlabel format(%9.2f)
    mlabsize(medium)         
    mlabposition(3 9 3 9)     
    mlabgap(*1.5)            
    mlabcolor(black)         
    
    graphregion(color(white) margin(small))
    bgcolor(white)
    plotregion(color(white) margin(medium))
    name(plot_direct_3, replace);
#delimit cr

// 4. Invalid Missing Responses
#delimit ;
coefplot 
    (ols_z_invalid_skips_c, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
        msymbol(circle) msize(medium) label("OLS"))
    (pds_z_invalid_skips_c, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
        msymbol(diamond) msize(medium) label("PDS-Lasso")),
    
    vertical
    xtitle("")
    xscale(range(0.5 2.5))
    xlabel(
        1 "T₁" 
        2 "T₂", 
        labsize(medium) angle(0)
    )
    
    yline(0, lpattern(dash) lcolor(gs10))
    yscale(range(-1.5 0))
    ylabel(-1.5(0.5)0, format(%9.2f) labsize(small))
    ytitle("Treatment Effect", size(medium))
    title("d) Invalid Missing Responses (z-score)", size(medsmall) color(black))
    
    legend(off)
    
    mlabel format(%9.2f)
    mlabsize(medium)          
    mlabposition(3 9 3 9)     
    mlabgap(*1.5)            
    mlabcolor(black)          
    
    graphregion(color(white) margin(small))
    bgcolor(white)
    plotregion(color(white) margin(medium))
    name(plot_direct_4, replace);
#delimit cr

// 5. Total Inconsistencies
#delimit ;
coefplot 
    (ols_z_inconsistency, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
        msymbol(circle) msize(medium) label("OLS"))
    (pds_z_inconsistency, keep(t1 t2) mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
        msymbol(diamond) msize(medium) label("PDS-Lasso")),
    
    vertical
    xtitle("")
    xscale(range(0.5 2.5))
    xlabel(
        1 "T₁" 
        2 "T₂", 
        labsize(medium) angle(0)
    )
    
    yline(0, lpattern(dash) lcolor(gs10))
    yscale(range(-0.3 0.5))
    ylabel(-0.3(0.2)0.5, format(%9.2f) labsize(small))
    ytitle("Treatment Effect", size(medium))
    title("e) Total Inconsistencies (z-score)", size(medsmall) color(black))

    legend(order(1 "OLS" 3 "PDS-Lasso") position(6) rows(1) size(medium))
        
    mlabel format(%9.2f)
    mlabsize(medium)         
    mlabposition(3 9 3 9)   
    mlabgap(*1.5)             
    mlabcolor(black)         
    
    graphregion(color(white) margin(small))
    bgcolor(white)
    plotregion(color(white) margin(medium))
    name(plot_direct_5, replace);
#delimit cr

#delimit ;
graph combine plot_direct_1 plot_direct_2 plot_direct_3 plot_direct_4 plot_direct_5,
    col(3)
    row(2)
    imargin(small)
    graphregion(color(white) margin(small))
    note("Note: Solid line/circle = OLS, Dashed line/diamond = PDS-Lasso", size(vsmall))
    ysize(6)
    xsize(12)
    scale(1)
    name(data_quality_direct, replace);
#delimit cr

graph save "$figures/Figure_A1.gph", replace
graph export "$figures/Figure_A1.png", replace 




/*==============================================================================
			Figure A2: Coefficient Plot of Impact on Nutrition Status
==============================================================================*/
use "$root\Data\child_level_anthro_vill.dta", clear


// Set up child-level controls
tab birth_order, gen(birth_order_)
tab primary_occu_mother, gen(mother_occu_)
tab primary_occu_father, gen(father_occu_)

tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"

global child_d_contr "child_female_ birth_order_* whether_employed_mother mother_occu_* father_occu_*"
global child_c_contr "age_in_months_ yrs_schooling_mother highest_class_father total_siblings total_hh_member q18"

/*==============================================================================
           HANDLING MISSING VALUES FOR PDS LASSO
==============================================================================*/

global miss_child_d_contr ""
global dmiss_child_d_contr ""

foreach var of varlist $child_d_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_d_contr "$miss_child_d_contr miss_`var'"
    global dmiss_child_d_contr "$dmiss_child_d_contr dmiss_`var'"
}

global miss_child_c_contr ""
global dmiss_child_c_contr ""

foreach var of varlist $child_c_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_c_contr "$miss_child_c_contr miss_`var'"
    global dmiss_child_c_contr "$dmiss_child_c_contr dmiss_`var'"
}

// Create final control variable global for child-level analysis
global child_controls "$str_contr $miss_child_d_contr $dmiss_child_d_contr $miss_child_c_contr $dmiss_child_c_contr"


eststo clear

// Panel A: OLS models
foreach var of varlist whether_stunted whether_wasted whether_underweight {
    areg `var' t1 t2 if child_age_ <=60, a(c_strata) vce(cluster vill_id)   
    eststo ols_`var'
}

// Panel B: PDS-Lasso models
foreach var of varlist whether_stunted whether_wasted whether_underweight {
    pdslasso `var' t1 t2 ($child_controls) if child_age_ <=60, cluster(vill_id) partial($str_contr)
    eststo pds_`var'
}


local outcomes "whether_stunted whether_wasted whether_underweight"
local titles `" "a) Whether Stunted" "b) Whether Wasted" "c) Whether Underweight" "'


local i = 1
foreach var of local outcomes {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "OLS" 3 "PDS-Lasso") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (ols_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
            msymbol(circle) msize(medium) label("OLS"))
        (pds_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
            msymbol(diamond) msize(medium) label("PDS-Lasso")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-0.15 0.1))
        ylabel(-0.15 (0.25)0.1, format(%9.2f) labsize(small))
        ytitle("Treatment Effect", size(small))
        title("`title'", size(small) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 9 3 9)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(anthro_plot`i', replace);
    #delimit cr
    
    local i = `i' + 1
}


#delimit ;
graph combine anthro_plot1 anthro_plot2 anthro_plot3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(small))
    note("Note: Solid line/circle = OLS, Dashed line/diamond = PDS-Lasso", size(vsmall))
    ysize(4)
    xsize(12)
    scale(1)
    name(anthro_combined, replace);
#delimit cr

graph save "$figures/Figure_A2.gph", replace
graph export "$figures/Figure_A2.png", replace width(2100) height(1400)








/*==============================================================================
			Figure A3: Coefficient Plot of Impact on Nutrition z-scores
==============================================================================*/
use "$root\Data\child_level_anthro_vill.dta", clear


// Set up child-level controls
tab birth_order, gen(birth_order_)
tab primary_occu_mother, gen(mother_occu_)
tab primary_occu_father, gen(father_occu_)

tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"

global child_d_contr "child_female_ birth_order_* whether_employed_mother mother_occu_* father_occu_*"
global child_c_contr "age_in_months_ yrs_schooling_mother highest_class_father total_siblings total_hh_member q18"


/*==============================================================================
           HANDLING MISSING VALUES FOR PDS LASSO
==============================================================================*/

global miss_child_d_contr ""
global dmiss_child_d_contr ""

foreach var of varlist $child_d_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_d_contr "$miss_child_d_contr miss_`var'"
    global dmiss_child_d_contr "$dmiss_child_d_contr dmiss_`var'"
}

global miss_child_c_contr ""
global dmiss_child_c_contr ""

foreach var of varlist $child_c_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_c_contr "$miss_child_c_contr miss_`var'"
    global dmiss_child_c_contr "$dmiss_child_c_contr dmiss_`var'"
}

// Create final control variable global for child-level analysis
global child_controls "$str_contr $miss_child_d_contr $dmiss_child_d_contr $miss_child_c_contr $dmiss_child_c_contr"




eststo clear

// Panel A: OLS models for z-scores
foreach var of varlist haz06 waz06 whz06 {
    areg `var' t1 t2 if child_age_ <=60 & (`var' >=-4 & `var' <=4), a(c_strata) vce(cluster vill_id)
    eststo ols_`var'
}

// Panel B: PDS-Lasso models for z-scores
foreach var of varlist haz06 waz06 whz06 {
    pdslasso `var' t1 t2 ($child_controls) if child_age_ <=60 & (`var' >=-4 & `var' <=4), cluster(vill_id) partial($str_contr)
    eststo pds_`var'
}



// Define properties
local outcomes "haz06 waz06 whz06"
local titles `" "a) Height-for-Age Z-score" "b) Weight-for-Age Z-score" "c) Weight-for-Height Z-score" "'

local i = 1
foreach var of local outcomes {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "OLS" 3 "PDS-Lasso") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    #delimit ;
    coefplot 
        (ols_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
            msymbol(circle) msize(small) label("OLS"))
        (pds_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
            msymbol(diamond) msize(small) label("PDS-Lasso")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(small) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-1.0 0.5))
        ylabel(-1.0(0.25)0.5, format(%9.2f) labsize(small))
        ytitle("Treatment Effect", size(vsmall))
        title("`title'", size(small) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 9 3 9)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(zscore_plot`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

#delimit ;
graph combine zscore_plot1 zscore_plot2 zscore_plot3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(small))
    note("Note: Solid line/circle = OLS, Dashed line/diamond = PDS-Lasso", size(vsmall))
    ysize(4)
    xsize(12)
    scale(1)
    name(child_anthropometric_zscores, replace);
#delimit cr

graph save "$figures/Figure_A3.gph", replace
graph export "$figures/Figure_A3.png", replace width(2100) height(1400)








/*==============================================================================
			Figure A4: Coefficient Plot of Impact on Child Consumption
==============================================================================*/
use "$root\Data\child_calorie.dta", clear

// Create education dummies if not already created
tab sec3_q11_, gen(education_)

// Define control variable groups
global child_cov_d "gender education_*"
global child_cov_c "child_age"
global hh_cov_d "father_illiterate father_pirimary_edu father_uppirimary_edu father_secondary_edu father_highsecondary_edu father_higher_edu mother_illiterate mother_pirimary_edu mother_uppirimary_edu mother_secondary_edu mother_highsecondary_edu mother_higher_edu ls_cas"
global hh_cov_c "hh_income mother_inc_cas"

// Create strata dummies
tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"

global all_covars "$child_cov_d $child_cov_c $hh_cov_d $hh_cov_c"

/*==============================================================================
           HANDLING MISSING VALUES FOR PDS LASSO 
==============================================================================*/

// First, handle missing values for discrete child controls
global miss_child_cov_d ""
global dmiss_child_cov_d ""

foreach var of varlist $child_cov_d {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_cov_d "$miss_child_cov_d miss_`var'"
    global dmiss_child_cov_d "$dmiss_child_cov_d dmiss_`var'"
}

// Then handle missing values for continuous child controls
global miss_child_cov_c ""
global dmiss_child_cov_c ""

foreach var of varlist $child_cov_c {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_child_cov_c "$miss_child_cov_c miss_`var'"
    global dmiss_child_cov_c "$dmiss_child_cov_c dmiss_`var'"
}

// Handle missing values for discrete household controls
global miss_hh_cov_d ""
global dmiss_hh_cov_d ""

foreach var of varlist $hh_cov_d {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_cov_d "$miss_hh_cov_d miss_`var'"
    global dmiss_hh_cov_d "$dmiss_hh_cov_d dmiss_`var'"
}

// Handle missing values for continuous household controls
global miss_hh_cov_c ""
global dmiss_hh_cov_c ""

foreach var of varlist $hh_cov_c {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_cov_c "$miss_hh_cov_c miss_`var'"
    global dmiss_hh_cov_c "$dmiss_hh_cov_c dmiss_`var'"
}

// Create consolidated control variable global for child consumption analysis
global control "$str_contr $miss_child_cov_d $dmiss_child_cov_d $miss_child_cov_c $dmiss_child_cov_c $miss_hh_cov_d $dmiss_hh_cov_d $miss_hh_cov_c $dmiss_hh_cov_c"


// Create shorter variable names, otherwise it is being truncated while stroring
gen calories = ch_calories_kcal_win_top20
gen protein = ch_protein_win_top20
gen fat = ch_fat_win_top20
gen carbs = ch_carbohydrates_win_top20

eststo clear

foreach var of varlist calories protein fat carbs {
    areg `var' t1 t2 child_age hh_income i.education_father i.education_mother, a(c_strata) vce(cluster vill_id)
    eststo ols_`var'
}

foreach var of varlist calories protein fat carbs {
    pdslasso `var' t1 t2 ($control), cluster(vill_id) partial($str_contr)
    eststo pds_`var'
}



local outcomes "calories protein fat carbs"
local titles `" "a) Daily Calorie (kcal)" "b) Protein (g)" "c) Fat (g)" "d) Carbohydrates (g)" "'
local names "Calorie Protein Fat Carbohydrates"

// Y-axis properties for each plot
local ymin1 = 0
local ymax1 = 150
local ystep1 = 30
local format1 = "%9.1f"

local ymin2 = 5
local ymax2 = 30
local ystep2 = 5
local format2 = "%9.1f"

local ymin3 = -2
local ymax3 = 10
local ystep3 = 2
local format3 = "%9.1f"

local ymin4 = 40
local ymax4 = 90
local ystep4 = 10
local format4 = "%9.1f"

local i = 1
foreach var of local outcomes {
    
    local title: word `i' of `titles'
    local name: word `i' of `names'
    local ymin = `ymin`i''
    local ymax = `ymax`i''
    local ystep = `ystep`i''
    local format = "`format`i''"
    
    // Legend only on last plot
    if `i' == 4 {
        local show_legend "legend(order(1 "OLS" 3 "PDS-Lasso") position(6) rows(1) size(small))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (ols_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
            msymbol(circle) msize(small) label("OLS"))
        (pds_`var', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
            msymbol(diamond) msize(small) label("PDS-Lasso")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(small) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(`ymin' `ymax'))
        ylabel(`ymin'(`ystep')`ymax', format(`format') labsize(small))
        ytitle("Treatment Effect", size(small))
        title("`title'", size(small) color(black))
        
        `show_legend'
        
        mlabel format(`format')
        mlabsize(small)
        mlabposition(3 9 3 9)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(`name', replace);
    #delimit cr
    
    local i = `i' + 1
}


#delimit ;
graph combine Calorie Protein Fat Carbohydrates,
    col(2)
    row(2)
    imargin(small)
    graphregion(color(white) margin(small))
    note("Note: Solid line/circle = OLS, Dashed line/diamond = PDS-Lasso", size(vsmall))
    ysize(8)
    xsize(12)
    scale(1)
    name(consumption_combined, replace);
#delimit cr

graph save "$figures/Figure_A4.gph", replace
graph export "$figures/Figure_A4.png", replace width(2100) height(1400)








/*==============================================================================
			Figure A5: Coefficient Plot of Impact on Psychological Well-being
==============================================================================*/
use "$root\Data\hh_dsy_cr_rec.dta", clear

// Religion
tab sec2_q1, gen(religion_)
// Caste
tab sec2_q2, gen(caste_)
// Household ownership type
tab sec2_q6, gen(house_own_type_)
// Other binary indicators
gen dbt_receipt = (sec2_q4==1)
gen improved_water = inlist(sec2_q18, 9, 10, 13, 14)
gen toilet_ownership = (sec2_q20==1)
gen improved_sanitation = inlist(sec2_q21, 2, 3, 4, 5) if sec2_q20 == 1
gen electricity = (sec2_q16==1)
gen internet = (sec2_q17==1)
gen clean_cooking = 0
replace clean_cooking = 1 if sec2_q24_1 == 1 | sec2_q24_2 == 1 | sec2_q24_3 == 1
gen child_mortality = (sec2_q23==1)
gen mother_working = (sec6_q3==1)

// Asset variables
egen asset_count = rowtotal(sec2_q12_1-sec2_q12_14 sec2_q12_96)
gen asset_proportion = asset_count / 15
egen asset_zscore = std(asset_proportion)
egen max_educ_hh = rowmax(sec3_q11_*)
gen hhsize = sec3_q1
gen num_rooms = sec2_q7
gen property_value = sec2_q8
gen dbt_amount = sec2_q5
egen total_health_exp = rowtotal(sec3_q18_*), missing
egen total_edu_exp = rowtotal(sec3_q17_*), missing
gen siblings_under12 = sec7_q2

// Define control variable groups
global hh_d_contr "religion_* caste_* house_own_type_* dbt_receipt electricity internet improved_water toilet_ownership clean_cooking child_mortality mother_working"
global hh_c_contr "hhsize max_educ_hh num_rooms property_value total_health_exp total_edu_exp asset_zscore siblings_under12"

tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"

/*==============================================================================
HANDLING MISSING VALUES FOR PDS LASSO, following David McKenzie's recommendations 
==============================================================================*/

global miss_hh_d_contr ""
global dmiss_hh_d_contr ""

foreach var of varlist $hh_d_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_d_contr "$miss_hh_d_contr miss_`var'"
    global dmiss_hh_d_contr "$dmiss_hh_d_contr dmiss_`var'"
}

global miss_hh_c_contr ""
global dmiss_hh_c_contr ""

foreach var of varlist $hh_c_contr {
    cap drop dmiss_`var'
    gen dmiss_`var' = missing(`var')
    
    cap drop miss_`var'
    gen miss_`var' = `var'
    replace miss_`var' = 0 if missing(`var')
    
    global miss_hh_c_contr "$miss_hh_c_contr miss_`var'"
    global dmiss_hh_c_contr "$dmiss_hh_c_contr dmiss_`var'"
}

// Create final control variable global
global controls "$str_contr $miss_hh_d_contr $dmiss_hh_d_contr $miss_hh_c_contr $dmiss_hh_c_contr"


zscore sec6_q47 sec6_q48 sec6_q49 sec6_q50 sec6_q51 sec6_q52 sec6_q53 sec6_q54 sec6_q55 sec6_q56

// Stress index 1 based on Cohen's stress scale (negative items)
gen stress1 = (z_sec6_q47+z_sec6_q48+z_sec6_q49+z_sec6_q52+z_sec6_q55+z_sec6_q56)/6

// Stress index 2 based on positive tone items
gen stress2 = (z_sec6_q50+ z_sec6_q51+ z_sec6_q53+ z_sec6_q54)/4

eststo clear

// Panel A: OLS models
foreach var of varlist stress1 stress2 var_wellbeing {
    areg `var' t1 t2 i.enum_id, a(c_strata) vce(cluster vill_id)
    eststo ols_`var'
}

// Panel B: PDS-Lasso models
foreach var of varlist stress1 stress2 var_wellbeing {
    pdslasso `var' t1 t2 ($controls), cluster(vill_id) partial($str_contr)
    eststo pds_`var'
}




// Define properties for each plot
local models "stress1 stress2 var_wellbeing"
local titles `" "a) Stress Index 1 (Negative items)" "b) Stress Index 2 (Positive items)" "c) Variance in Psychological Wellbeing" "'

// Y-axis properties
local ymin1 = -2.4
local ymax1 = 0.2
local ystep1 = 0.5

local ymin2 = -0.4
local ymax2 = 3.0
local ystep2 = 1.0

local ymin3 = -0.4
local ymax3 = 0.8
local ystep3 = 0.4

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    local ymin = `ymin`i''
    local ymax = `ymax`i''
    local ystep = `ystep`i''
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "OLS" 3 "PDS-Lasso") position(6) rows(1) size(small))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (ols_`model', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(medthick)) 
            msymbol(circle) msize(small) label("OLS"))
        (pds_`model', keep(t1 t2) mcolor(black) 
            ciopts(lcolor(black) recast(rcap) lwidth(thick) lpattern(dash)) 
            msymbol(diamond) msize(small) label("PDS-Lasso")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(small) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(`ymin' `ymax'))
        ylabel(`ymin'(`ystep')`ymax', format(%9.2f) labsize(small))
        ytitle("Treatment Effect", size(small))
        title("`title'", size(small) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(small)
        mlabposition(3 9 3 9)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(small))
        name(stress_plot`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine stress_plot1 stress_plot2 stress_plot3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(vsmall))
    note("Note: Solid line/circle = OLS, Dashed line/diamond = PDS-Lasso", size(vsmall))
    ysize(4)
    xsize(12)
    scale(1)
    name(stress_combined, replace);
#delimit cr

graph save "$figures/Figure_A5.gph", replace
graph export "$figures/Figure_A5.png", replace width(1600) height(800)






/*==============================================================================
			Appendix 10 (Figure A6-Figure A8) Robustness Check: Benford's Law
==============================================================================*/

/*==========================================================================================
Figure A6: Robustness Check using Benford's Law: Coefficient Plot of Child Nutrition Status
==========================================================================================*/
use "$root\Data\benford_anthro.dta", clear

// Run all regressions
eststo clear
eststo m1: reg ndist_height t1 t2 q4 q15 q19, robust
eststo m2: reg ndist_weight t1 t2 q4 q15 q19, robust
eststo m3: reg ndist_child_age_ t1 t2 q4 q15 q19, robust
eststo m4: reg ndist_haz06 t1 t2 q4 q15 q19, robust
eststo m5: reg ndist_waz06 t1 t2 q4 q15 q19, robust
eststo m6: reg ndist_whz06 t1 t2 q4 q15 q19, robust

// Define properties
local models "m1 m2 m3 m4 m5 m6"
local titles `" "a) Height" "b) Weight" "c) Child Age" "d) HAZ" "e) WAZ" "f) WHZ" "'

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 6 {
        local show_legend "legend(order(1 "T₁: Retesting" 3 "T₂: Agentic Feedback") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (`model', keep(t1) rename(t1=1)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(solid))
            msymbol(circle) msize(medium) label("T₁: Retesting"))
        (`model', keep(t2) rename(t2=2)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(dash))
            msymbol(diamond) msize(medium) label("T₂: Agentic Feedback")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-1.0 2.0))
        ylabel(-1.0(0.5)2.0, format(%9.1f) labsize(small))
        ytitle("Deviation from Benford's Distribution", size(vsmall))
        title("`title'", size(medsmall) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 3)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(anthro_plot_`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine anthro_plot_1 anthro_plot_2 anthro_plot_3 
              anthro_plot_4 anthro_plot_5 anthro_plot_6,
    col(3)
    row(2)
    imargin(small)
    graphregion(color(white) margin(small))
    ysize(6)
    xsize(14)
    scale(1)
    name(benford_anthro_combined, replace);
#delimit cr

graph save "$figures/Figure_A6.gph", replace
graph export "$figures/Figure_A6.png", replace width(4200)










/*==========================================================================================
Figure A7: Robustness Check using Benford's Law: Coefficient Plot of Child Consumption
==========================================================================================*/
use "$root\Data\benford_con.dta", clear


// Run all regressions
eststo clear
eststo m1: reg ndist_protein t1 t2 q4 q15 q19, robust
eststo m2: reg ndist_carbs t1 t2 q4 q15 q19, robust
eststo m3: reg ndist_cal2 t1 t2 q4 q15 q19, robust
eststo m4: reg ndist_fat t1 t2 q4 q15 q19, robust

// Define properties
local models "m1 m2 m3 m4"
local titles `" "a) Protein" "b) Carbohydrates" "c) Calories" "d) Fat" "'

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 4 {
        local show_legend "legend(order(1 "T₁: Retesting" 3 "T₂: Agentic Feedback") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (`model', keep(t1) rename(t1=1)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(solid))
            msymbol(circle) msize(medium) label("T₁: Retesting"))
        (`model', keep(t2) rename(t2=2)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(dash))
            msymbol(diamond) msize(medium) label("T₂: Agentic Feedback")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-2.2 0.5))
        ylabel(-2.0(0.5)0.5, format(%9.1f) labsize(small))
        ytitle("Deviation from Benford's Distribution", size(vsmall))
        title("`title'", size(medsmall) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 3)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(benford_plot_`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine benford_plot_1 benford_plot_2 benford_plot_3 benford_plot_4,
    col(2)
    row(2)
    imargin(small)
    graphregion(color(white) margin(small))
    ysize(8)
    xsize(16)
    scale(1)
    name(benford_nutrition_combined, replace);
#delimit cr

graph save "$figures/Figure_A7.gph", replace
graph export "$figures/Figure_A7.png", replace width(4800)


/*==========================================================================================
Figure A8: Robustness Check using Benford's Law: Coefficient Plot of Stress Score
==========================================================================================*/


use "$root\Data\benford_stress.dta", clear

// Run all regressions
eststo clear
eststo m1: reg ndist_stress_score1 t1 t2 q4 q15 q19, robust
eststo m2: reg ndist_stress_score2 t1 t2 q4 q15 q19, robust 
eststo m3: reg ndist_stress_score_comp t1 t2 q4 q15 q19, robust

// Define properties
local models "m1 m2 m3"
local titles `" "a) Stress Score 1" "b) Stress Score 2" "c) Composite Score" "'

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "T₁: Retesting" 3 "T₂: Agentic Feedback") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (`model', keep(t1) rename(t1=1)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(solid))
            msymbol(circle) msize(medium) label("T₁: Retesting"))
        (`model', keep(t2) rename(t2=2)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(dash))
            msymbol(diamond) msize(medium) label("T₂: Agentic Feedback")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-2.0 1.0))
        ylabel(-2.0(0.5)1.0, format(%9.1f) labsize(small))
        ytitle("Deviation from Benford's Distribution", size(small))
        title("`title'", size(medsmall) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 3)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(stress_plot_`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine stress_plot_1 stress_plot_2 stress_plot_3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(small))
    ysize(4)
    xsize(12)
    scale(1)
    name(benford_stress_combined, replace);
#delimit cr

graph save "$figures/Figure_A8.gph", replace
graph export "$figures/Figure_A8.png", replace width(1600) height(800)







/*==============================================================================
		Appendix 11: Robustness Check: Enumerator Effect (Figure A9-Figure A10)
==============================================================================*/
/*==============================================================================
			Figure A9: Enumerator Effect and Child Nutrition Indicators
================================================================================*/
use "$root\Data\child_calorie.dta", clear
sort vill_id
egen vill_id_group = group(vill_id)
tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"
global nutrition "ch_calories_kcal_win_top20 ch_protein_win_top20 ch_fat_win_top20 ch_carbohydrates_win_top20"
local i = 1
foreach var of global nutrition {
    // Treatment 1 (Retesting)
    reg `var' i.enum_id if t1 == 1, vce(cluster vill_id)
    gen rsq_t1_v`i' = e(r2)
    
    // Treatment 2 (Agentic Feedback)
    reg `var' i.enum_id if t2 == 1, vce(cluster vill_id)
    gen rsq_t2_v`i' = e(r2)
    
    // Control group (Restudying)
    reg `var' i.enum_id if t1 == 0 & t2 == 0, vce(cluster vill_id)
    gen rsq_control_v`i' = e(r2)
    
    local ++i
}


preserve
keep rsq_*
keep in 1
gen id = 1
reshape long rsq_t1_v rsq_t2_v rsq_control_v, i(id) j(varnum)
rename rsq_t1_v rsq_t1_val
rename rsq_t2_v rsq_t2_val  
rename rsq_control_v rsq_control_val
gen variable = ""
replace variable = "ch_calories_kcal_win_top20" if varnum == 1
replace variable = "ch_protein_win_top20" if varnum == 2
replace variable = "ch_fat_win_top20" if varnum == 3
replace variable = "ch_carbohydrates_win_top20" if varnum == 4
gen order = varnum
gen short_label = ""
replace short_label = "Daily Calorie" if varnum == 1
replace short_label = "Protein" if varnum == 2
replace short_label = "Fat" if varnum == 3
replace short_label = "Carbohydrates" if varnum == 4
list variable short_label rsq_t1_val rsq_t2_val rsq_control_val order
sort order



#delimit ;
twoway 
    // Retesting (T1) line and markers - using black color with solid line
    (connected rsq_t1_val order,
        msymbol(circle) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(solid))
        
    // Agentic Feedback (T2) line and markers - using black color with dashed line
    (connected rsq_t2_val order,
        msymbol(diamond) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(dash))
        
    // Restudying (Control) line and markers - using black color with dotted line
    (connected rsq_control_val order,
        msymbol(square) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(dot)),
        
    ytitle("R-Squared", size(small)) 
    xtitle("Child Nutrition Indicators", size(small))
    ylabel(0.05(0.03)0.3, angle(horizontal) labsize(small) grid glcolor(gs14) glwidth(vvthin)) 
    yscale(range(0.05 0.3))
    
    xlabel(1 `""Daily" "Calorie""' 
           2 `"Protein"' 
           3 `"Fat"' 
           4 `"Carbohydrates"', 
           labsize(small) labgap(2))
    
    legend(order(
        1 "Retesting" 
        2 "Agentic Feedback" 
        3 "Restudying"
    ) rows(1) position(6) size(small))
        
    graphregion(color(white) margin(large)) 
    plotregion(margin(medium))
    bgcolor(white)
    xsize(6.5) ysize(4.5)
    name(child_nutrition_rsq_plot, replace);
#delimit cr


graph save "$figures/Figure_A9.gph", replace
graph export "$figures/Figure_A9.png", replace
restore







/*==============================================================================
			Figure A10: Enumerator Effect and Stress Items and Index
================================================================================*/
use "$root\Data\hh_dsy_cr_rec.dta", clear
sort vill_id
egen vill_id_group = group(vill_id)
tab c_strata, gen(c_strata_)
global str_contr "c_strata_*"

zscore sec6_q50 sec6_q51 sec6_q53 sec6_q54 


gen stress_index2 = (z_sec6_q50 + z_sec6_q51 + z_sec6_q53 + z_sec6_q54) / 4  	 // Positive experiences/coping items
label variable z_sec6_q50 "Confident handling personal problems"
label variable z_sec6_q51 "Things going your way"
label variable z_sec6_q53 "Able to control irritations"
label variable z_sec6_q54 "Felt on top of things"
label variable stress_index2 "Stress Index (Positive Items)"



// Define positive coping variables and their index
global positive_coping "z_sec6_q50 z_sec6_q51 z_sec6_q53 z_sec6_q54 stress_index2"

preserve

foreach var of global positive_coping {
    // Treatment 1 (Retesting)
    reg `var' i.enum_id if t1 == 1, vce(cluster vill_id)
    gen rsq_t1_`var' = e(r2)
    
    // Treatment 2 (Agentic Feedback)
    reg `var' i.enum_id if t2 == 1, vce(cluster vill_id)
    gen rsq_t2_`var' = e(r2)
    
    // Control group (Restudying)
    reg `var' i.enum_id if t1 == 0 & t2 == 0, vce(cluster vill_id)
    gen rsq_control_`var' = e(r2)
}

keep rsq_*
keep in 1
gen id = 1

reshape long rsq_t1_ rsq_t2_ rsq_control_, i(id) j(variable) string
rename rsq_t1_ rsq_t1_val
rename rsq_t2_ rsq_t2_val  
rename rsq_control_ rsq_control_val

gen order = .
replace order = 1 if variable == "z_sec6_q50"       // "Confident handling personal problems"
replace order = 2 if variable == "z_sec6_q51"       // "Things going your way"
replace order = 3 if variable == "z_sec6_q53"       // "Able to control irritations"
replace order = 4 if variable == "z_sec6_q54"       // "Felt on top of things"
replace order = 5 if variable == "stress_index2"    // "Stress Index (Positive Items)"

gen short_label = ""
replace short_label = "Confident handling personal problems" if variable == "z_sec6_q50"
replace short_label = "Things going your way" if variable == "z_sec6_q51"
replace short_label = "Able to control irritations" if variable == "z_sec6_q53"
replace short_label = "Felt on top of things" if variable == "z_sec6_q54"
replace short_label = "Stress Index (Positive Items)" if variable == "stress_index2"

list variable short_label rsq_t1_val rsq_t2_val rsq_control_val order
sort order

#delimit ;
twoway 
    // Retesting (T1) line and markers - using black color with solid line
    (connected rsq_t1_val order,
        msymbol(circle) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(solid))
        
    // Agentic Feedback (T2) line and markers - using black color with dashed line
    (connected rsq_t2_val order,
        msymbol(diamond) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(dash))
        
    // Restudying (Control) line and markers - using black color with dotted line
    (connected rsq_control_val order,
        msymbol(square) mcolor(black) msize(small) 
        lcolor(black) lwidth(thin) lpattern(dot)),
        
    ytitle("R-Squared", size(small)) 
    xtitle("Positive Coping Items and Index", size(small))
    ylabel(0.1(0.05)0.6, angle(horizontal) labsize(small) grid glcolor(gs14) glwidth(vvthin)) 
    yscale(range(0.1 0.6))
    
    xlabel(1 `""Confident handling" "personal problems""' 
           2 `""Things going" "your way""' 
           3 `""Able to control" "irritations""' 
           4 `""Felt on top" "of things""'
           5 `""Stress Index 2" "(Positive Items)""', 
           labsize(small) labgap(2))
    
    legend(order(
        1 "Retesting" 
        2 "Agentic Feedback" 
        3 "Restudying"
    ) rows(1) position(6) size(small))
        
    graphregion(color(white) margin(large)) 
    plotregion(margin(medium))
    bgcolor(white)
    xsize(7) ysize(4.5)
    name(positive_coping_rsq_plot, replace);
#delimit cr

graph export "$figures/Figure_A10.png", replace
restore





/*==============================================================================
    Appendix 13: Mechanism - Self-Confidence, Motivation, Self-Responsibility
==============================================================================*/
use "$root\Data\enum_train.dta", clear

eststo clear

// Run regressions
eststo m1: areg list_score1 t1 t2 if list_treat==1, a(c_strata) robust
eststo m2: areg list_score2 t1 t2 if list_treat==1, a(c_strata) robust
eststo m3: areg list_score3 t1 t2 if list_treat==1, a(c_strata) robust

// Define properties
local models "m1 m2 m3"
local titles `" "a) Self Responsibility" "b) Self Confidence" "c) Motivation" "'

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "T₁: Retesting" 3 "T₂: Agentic Feedback") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (`model', keep(t1) rename(t1=1)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(solid))
            msymbol(circle) msize(medium) label("T₁: Retesting"))
        (`model', keep(t2) rename(t2=2)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(dash))
            msymbol(diamond) msize(medium) label("T₂: Agentic Feedback")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(-3 1))
        ylabel(-3(1)1, format(%9.2f) labsize(small))
        ytitle("Treatment Effect", size(small))
        title("`title'", size(medsmall) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 3)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(list_plot_`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine list_plot_1 list_plot_2 list_plot_3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(small))
    ysize(4)
    xsize(12)
    scale(1)
    name(list_experiment_combined, replace);
#delimit cr

graph save "$figures/Appendix_13.gph", replace
graph export "$figures/Appendix_13.png", replace width(1600) height(800)








/*==============================================================================
Appendix 14: Mechanism: My Enquiries Reveal My Methods
================================================================================*/

use "$root\Data\enum_train.dta", clear

eststo clear
eststo m1: areg prop_q1 t1 t2, a(c_strata) robust
eststo m2: areg prop_q2 t1 t2, a(c_strata) robust
eststo m3: areg prop_q3 t1 t2, a(c_strata) robust

// Define properties
local models "m1 m2 m3"
local titles `" "a) Anticipatory Questions" "b) Inferential Questions" "c) Confirmatory Questions" "'

// Define y-axis properties (different for each plot)
local ymin1 = -0.8
local ymax1 = 0.2
local ystep1 = 0.2

local ymin2 = -0.2
local ymax2 = 0.6
local ystep2 = 0.2

local ymin3 = -0.2
local ymax3 = 0.6
local ystep3 = 0.2

local i = 1
foreach model of local models {
    
    local title: word `i' of `titles'
    local ymin = `ymin`i''
    local ymax = `ymax`i''
    local ystep = `ystep`i''
    
    // Legend only on last plot
    if `i' == 3 {
        local show_legend "legend(order(1 "T₁: Retesting" 3 "T₂: Agentic Feedback") position(6) rows(1) size(medium))"
    }
    else {
        local show_legend "legend(off)"
    }
    
    // Create plot
    #delimit ;
    coefplot 
        (`model', keep(t1) rename(t1=1)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(solid))
            msymbol(circle) msize(medium) label("T₁: Retesting"))
        (`model', keep(t2) rename(t2=2)
            mcolor(black) ciopts(lcolor(black) recast(rcap) lwidth(medthick) lpattern(dash))
            msymbol(diamond) msize(medium) label("T₂: Agentic Feedback")),
        
        vertical
        xtitle("")
        xscale(range(0.5 2.5))
        xlabel(1 "T₁" 2 "T₂", labsize(medium) angle(0))
        
        yline(0, lpattern(dash) lcolor(gs10))
        yscale(range(`ymin' `ymax'))
        ylabel(`ymin'(`ystep')`ymax', format(%9.2f) labsize(small))
        ytitle("Treatment Effect", size(small))
        title("`title'", size(medsmall) color(black))
        
        `show_legend'
        
        mlabel format(%9.2f)
        mlabsize(medium)
        mlabposition(3 3)
        mlabgap(*1.5)
        mlabcolor(black)
        
        graphregion(color(white) margin(small))
        bgcolor(white)
        plotregion(color(white) margin(medium))
        name(call_plot_`i', replace);
    #delimit cr
    
    local i = `i' + 1
}

// Combine plots
#delimit ;
graph combine call_plot_1 call_plot_2 call_plot_3,
    col(2)
    row(1)
    imargin(small)
    graphregion(color(white) margin(small))
    ysize(4)
    xsize(12)
    scale(1)
    name(call_logs_combined, replace);
#delimit cr

graph save "$figures/Appendix_14.gph", replace
graph export "$figures/Appendix_14.png", replace width(1600) height(800)




